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Abstract. An algorithm for the computer simulation of large space structures 
under active control is considered. Linked lists are used in a matrix data structure 
to implement the trapezoidal rule on the system differential equations. The use of the 
trapezoidal rule ensures that the numerical stability is equivalent to the system stability, 
which is essential for this type of simulation. The sparsity of the system matrices is 
exploited by the linked lists, and the algorithm efficiently steps through the lists in an 
orderly fashion. Results of simulations on a NASA large space structure experiment are 
reported. 
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1. Introduction. Future large space structures such as the space crane and the 
aerobrake will possess high flexibility and low damping. The suppression of vibrations is 
an important problem for this type of structure, since slewing maneuvers or disturbances 
can cause large amplitude vibrations over long time intervals. 

Vibration suppression can be accomplished via active feedback control. However, 
it is possible that an unstable controller design would damage the structure. Com- 
puter simulations are therefore desirable for evaluating controller performance and for 
detecting instability. 

Various approaches have been taken to simulate controller - structure interaction 
(CSI), such as those described in [1,2, 3, 4] and [5,p.3]. These approaches have accuracy 
and / or numerical stability limitations inherent in them. The method to be described 
in this paper overcomes these limitations. For instance, accuracy is maximized through 
the use of the finite element model of the large space structure, instead of a reduced 
order model. In addition, the numerical stability of the simulation is made equiva- 
lent to the stability of the physical system. This has been difficult to implement in 
the past, because of the differences between large space structure models and control 
system models. Finite element models of large space structures contain large, sparse, 
and symmetric matrices; compensator models tend to be small, dense, and unsymmet- 
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ric. An implicit integration scheme such as the trapezoidal rule has the property that 
the numerical stability of the simulation is equivalent to the stability of the physical 
system, but it also combines the matrices of the structure and the compensator; this 
destroys sparsity and symmetry. Computer memory requirements can thereby become 
unpractically large, if finite element data structures are used (banded matrix, skyline 
matrix, etc.). Therefore, a linked list data structure is developed below which recovers 
the sparsity of the computer simulation. Linked list algebra and factorization are also 
developed, so that an implicit integration scheme can be implemented. 

2. Closed Loop System Equations. Consider the following linear finite element 
model of a multi - input, multi - output structure: 


Mq + Cq + Kq = Bu (1) 

/ q \ 

y = C y I q (2) 

\q / 

where the displacement q e $ff nxl , the velocity q e 5J nxl , the acceleration q e 5R nXl , the 
input force u e 5ft mxl , the output vector y e 3? 3nxl , and M, C, K, B, C y are the mass, 
damping, stiffness, input, and output matrices, respectively. 

The mass matrix is assumed to be positive definite; the damping and stiffness ma- 
trices are assumed to be positive semidefinite. An additional n differential equations 
are associated with equation (1), because the velocities are the derivatives of the dis- 
placements. 

A linear compensator is assumed, with dynamics as follows: 

(:)-(;) <*> 

where the compensator state vector x e S? rxl , and the matrix L e §fj( r +m)x(r+3n) 

The differential equations above can be appended together to form one set of differential 
equations: 

z = Vz (4) 

If we have the state z at the Nth time step (time t), then we can obtain the state 
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( 5 ) 


at the N+l time step (time t+h) by using the trapezoidal rule: 

z N+i = Z N 2 ^ n + x ^ *n) 

A linear system is stable if and only if its eigenvalues are in the left half complex 
plane. It is desirable that this system stability region coincide with the numerical stabil- 
ity region. Unexpected damage would result if the computer simulation of a controlled 
large space structure were to show stability, when in fact the controller were to destabi- 
lize the structure. The numerical stability region of the trapezoidal rule does coincide 
with the left half complex plane [6,7]; therefore this algorithm is a logical choice for the 
simulation of controller — structure interaction. 

3. Sparse Matrix Storage. Consider the following simplified control problem, 
which illustrates the computer storage difficulties associated with the simulation of 
controller - structure interaction: 

If the structure is in a steady state condition, equation (1) simplifies to: 


Kq = Bu 


( 6 ) 


As an example, consider the case where the structure has a tridiagonal stiffness matrix 
(the ‘x’ entries represent nonzero elements): 
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The nonzero elements of this matrix exhibit a certain pattern, which makes it possible 
to store the three diagonals of the matrix as three arrays. Let us place a force actuator 
at the first degree of freedom: 


B = 


/ 1 \ 
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VO } 


( 8 ) 
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Now place a displacement sensor at the last degree of freedom: 



y — Cq 1 ! 



( 9 ) 

C q - 

(0 0 0 0 

i) 

(10) 

and establish output feedback from the displacement sensor to the actuator: 
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II 

3 



(11) 

Then the closed loop system equation becomes: 





Kq = 0 



(12) 


f x x 0 0 





a; x x 0 

0 



K = 

0 x x x 
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(13) 
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The control has affected the sparsity pattern, and it is not clear if renumbering the 
degrees of freedom would help considerably. When implementing the trapezoidal rule on 
equations (1-3), the same situation is encountered. If all elements of such matrices were 
stored, then computer memory would be wasted on the storage of zero-valued matrix 
elements. This would be a serious problem for large space structure type problems. 
Therefore, a linked list data structure is developed below which stores only the nonzero 
elements. 

4. Linked List Matrices. Figure 1 displays the data structure for a linked list 
matrix. The leftmost array in the figure contains the number of nonzero columns in each 
row. Adjacent to this array is another array which points into linked lists for each row. 
Each record in the linked list contains a field for a floating point matrix data element, 
and a field for the column number. 

As it will be shown later, addition and multiplication of matrices can be done if 
the linked lists are traversed in one direction. However, factorization of matrices will 
require the deletion of matrix elements. A matrix element deletion requires that the 
two surrounding elements be connected together; knowledge of the locations of the two 
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(Arrays) 


Figure 1: Linked List Matrix 
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surrounding elements is needed. This information can be quickly obtained if the linked 
list can be traversed in both directions. To establish this “double linking”, there are 
two pointers in each linked list record which point to the two surrounding elements. If 
there is only one element in a linked list, then the two pointers point at that element. 

Although it is possible to store and recall elements in any order (random access), it 
is more efficient to store and recall matrices by rows. The next section demonstrates 
that such row access can be used to implement the algebra needed for a controller - 
structure interaction simulation. 

5. Sparse Matrix Algebra 

5.1 Sparse Matrix Addition and Multiplication 

The addition of matrices is necessary for the assembly of the finite element model 
from the element mass and stiffness matrices. Matrix addition, multiplication, and 
factorization is also required for the implementation of the trapezoidal rule (see [5] 
for details). The addition of two linked list matrices can be accomplished by stepping 
through the linked lists of the first matrix by rows. At the same time, the corresponding 
element in the second matrix is recalled. The two elements from the two matrices are 
added together and stored in the second matrix. Thus the final result appears in the 
second matrix. 

The multiplication of two matrices is often computed by using inner products: 

Given A e 9? m xn , Be9t nxp 

n 

C = A*B, Ceft mxp , C ij =^A ik B kj (14) 

k=l 

In order to perform the multiplication efficiently by stepping through the linked lists 
in order, it is necessary to rearrange formula (14) into a form which resembles an outer 
product: 

m n 

c = EE a / ( 15 ) 

i=l j=l 

where /3 tJ is the zero matrix of dimension (m X p ) with the ith row replaced by the 
jth row of the B matrix. Note that the A matrix is stepped through once, and that the 
B matrix is stepped through at most m times. This multiplication procedure can be 
illustrated for two dimensional matrices: 
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( + 0-12621 

\021611 + <122621 


<111612 + <112622 
021612 + 022622 



5.2 LU factorization of a Sparse Matrix 

It is known that a matrix A can be factorized into a product of a lower triangular 
matrix L and an upper triangular matrix U. If the matrix equation Ax = b must be 
solved repeatedly for different b vectors, then the LU factorization leads to computa- 
tional efficiency [8]. This situation exists during the simulation of controller - structure 
interaction (see [5]). A linked list matrix is to be LU factorized, and therefore the usual 
procedure for LU factorization needs to be modified. The new procedure is given as 
follows: 

a) Form the transpose of the matrix A (A T ). 

b) For all of the rows which have not been selected as a pivot row: 

1. Select the sparsest row as the pivot row. 

This idea was used in [9]. However, the procedure described in [9] assumed that 
the nonzero elements of the sparse matrix to be factorized are packed into arrays. 
Deletions of matrix elements are necessary in the factorization, which is cumber- 
some if the data is stored in arrays. With linked lists, deletions are simple in that 
pointers in the matrix elements surrounding the deleted element are redirected 
at each other. The insertion of matrix elements is also relatively simple. This is 
important for the assembly of the finite element model via the addition of element 
mass and stiffness matrices into the global mass and stiffness matrices [10]. Thus 
a single data structure (the linked list matrix) can be used for the assembly of 
the finite element model and for the implementation of the trapezoidal rule; this 
avoids the need for conversions between data structures. 


2. Find the element with the maximum magnitude within the pivot row. This ele- 
ment will be referred to as the pivot element, and the column of this element will 
be referred to as the pivot column. Later below, division by the pivot element 
will be performed, which is why the maximum magnitude element is chosen as 
the pivot element. If the matrix is nonsingular, a nonzero element will be found 
in the pivot row. 

3. Delete the pivot row from A T . 

4. Save a copy of the pivot element, and delete it from A. 

5. The nonzero elements in the pivot column are referred to as “target elements”. 
These target elements appear in a row in A T , allowing for quick access. For each 
target element, 

i) multiplier = - target element / pivot element 

ii) Store the multiplier in the matrix L, and delete the target element from A 

and A t 

iii) Multiply the pivot row by the multiplier and add it to the target row of A 
and A t Since we chose the sparsest row as the pivot row, we retain as much 
sparsity as possible. 

(End of factorization: The matrix U now appears in place of matrix A) 

The matrix A is stored by rows. Since access by columns is needed above, A A is 
stored. It might appear that this would double the storage requirements. However, 
A is sparse. The LU factorization process creates “fill in” (some of the zero elements 
become nonzero) within the U matrix. After each step of the factorization process, 
another column of A T is no longer needed. The storage for this column goes towards 
the storage of the “fill in” elements. 

After the factorization has been completed, the linked list matrix A is left with the 
upper triangular part of the factorization, and linked list matrix L is the lower triangular 
part. 

5.3 LDL t Factorization of a Sparse Matrix 

A symmetric linked list matrix equation needs to be solved at every time step in 
the simulation of controller - structure interaction (see [5]). Another problem where 
a symmetric linked list matrix equation has to be repeatedly solved occurs during the 
computation of the lowest frequency modes of a structure (see [5]). The LDL T factoriza- 
tion can be employed to efficiently solve these matrix equations. L is a lower triangular 
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matrix with ones on the diagonal, and D is a diagonal matrix. The LDL T factorization 
is described and analyzed in Golub and Van Loan [8]. The algorithm is listed below in 
the style which Golub and Van Loan use in their text. The notation l:j signifies all of 
the integers from 1 to j. 

LDL t Algorithm: If Ae3ft( nxn ) is symmetric then the following algorithm computes a 
unit lower triangular matrix L and a diagonal matrix D so that A = LDL^. It is 
assumed that only the lower half of the matrix A is stored, because of the symmetry. 
The matrix A is overwritten with the matrices L and D by this algorithm. 


for j = 1 : 7i 


for z = 1 : (j — 1) 
v(i) = A(j y i)A(i y i ) 

end 

v(j) = A{j,j)~ A(j,(l : j - l))t»(l : 
= v{j) 

A({j + 1) : n,j) = 

(A((j + 1) : n,j) - A((j + 1) : n, 1 


(3 ~ 1)) 


:(,--l))*(l:(j -!)))/*(*) 


The LDL t algorithm can be modified to handle sparse symmetric matrices of the linked 
list storage type: 

The formation of the v vector on lines (2 - 5) is done by stepping through linked 
lists, instead of looping over the entire range from 1 to (j - 1). Because of the sparsity 
of the matrix, the v vector is also sparse. 

The algorithm yields one column of the L matrix at a time, as shown in lines (7 - 
8) of the algorithm. In line (8), a vector is formed by multiplying a submatrix by the 
v vector (A((j + 1) : n, 1 : (j — l))v(l : (j - 1))). A small example will be useful for 
illustrating the sparsity: 


v r = Mvj 


(16) 
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The symbol u x n in the above equation signifies a nonzero element which is in a 
linked list. There is no need to perform the computations for rows 3 and 4 of the result 
vector v r , because all of the terms for those rows are zero. Since the linked list storage 
scheme is being used, the zero part of v; does not actually exist: 


v r = Mv] 
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(18) 


(19) 


The matrix M is stored by rows. If the transpose of M is stored, then it is efficient 
to traverse the columns of M corresponding to the nonzero rows in v t . The union of 
the nonzero rows in those columns forms the set of nonzero rows in the result v r . We 
will refer to this set of rows which need to be handled as the set |J. 

It might appear that the storage of A T would double the storage requirements. 
However, A is sparse. The LDL T factorization process creates “fill in” within the 
matrix. However, after each step of the factorization process, another column of A A is 
no longer needed. The storage for this column goes towards the storage of the “fill in” 
elements. 

The new algorithm is given as follows: 

a) Form the transpose of the matrix A. 

b) For each row p of A: 

1. Delete the row from the A . 

2. Determine the set U of rows which need to handled. 
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3. Recall the diagonal element A pp for that row, and delete it from A. 

4. Compute the v vector described above. 

5. For all the rows q in the set |J which need to be handled: 

i) Recall the section of the row q of A up to column p. 

ii) Compute the dot product of that section with the v vector. 

iii) Subtract this dot product from A qp , and divide by v p . Store this result in 

Aqp and in Apq. 

(End of Procedure) 

6. The Mini-Mast Truss. Mini-Mast is an active control experiment being 
maintained at the NASA Langley Research Center [11]. A linear finite element model 
having 714 degrees of freedom was developed for this truss. Two of the accelerometers 
were used as sensors, and all three of the torque wheels were used for control. 

The Rayleigh damping coefficients were tuned to provide 2 percent damping in the 
first two modes and 1 percent damping in the next three modes. These first five modes 
and the dynamics of the torque wheels were combined to form a reduced order model 
of the structure, and linear quadratic regulator theory was used to design a controller. 
The total number of states in the compensator state vector is 16. Figures 2 and 3 show 
the open loop response and the closed loop response, respectively. In both cases, a one 
second triangular pulse was applied at one of the torque wheels. The simulations were 
done on a Sun-4 workstation; total memory requirements were less than 4 megabytes. 
About 6 minutes of computer time was used to produce the 280 time steps shown in 
the simulation. 

7. Conclusion. An alternative approach towards the simulation of controller - 
structure interaction (CSI) has been described in this paper. Linked lists were used 
to implement the trapezoidal rule, which enforces an equivalence between numerical 
stability and system stability. This characteristic is essential for CSI analysis, and has 
not been demonstrated by previous CSI simulation methods. 

Matrix storage has been implemented with linked lists, which required the develop- 
ment of linked list matrix algebra for the implementation of the trapezoidal rule. Thus 
methods for linked list matrix addition, multiplication, LU factorization, and LDL T 
factorization were developed. The linked lists are stepped through in order in these 
methods, which minimizes the number of computations required. Memory requirements 
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Figure 2: Mini-Mast Open Loop Response 
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Figure 3: Mini-Mast Closed Loop Response 
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are also minimized because storage is dedicated only to nonzero matrix elements. 

The feasability of this approach was demonstrated on a computer workstation for 
a large space structure experiment (the Mini-Mast truss). These linked list matrix 
methods show that it is possible to simulate the control of some large space structures 
without the use of a supercomputer. In the case of a supercomputer, it is not certain how 
effective linked list methods would be. A linked list is not in the form of a vector, which 
suggests that methods based on it would not take advantage of the special capabilities of 
vector processing machines. In the case of parallel processing supercomputers, research 
needs to be done to determine the effectiveness of these methods on such machines. 
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